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Abstract 

CN 

5_! Laminated glass structures are formed by stiff layers of glass connected with a compliant 

plastic interlayer. Due to their slenderness and heterogeneity, they exhibit a complex me- 
chanical response that is difficult to capture by single-layer models even in the elastic range. 
The purpose of this paper is to introduce an efficient and reliable finite element approach to 
the simulation of the immediate response of laminated glass beams. It proceeds from a refined 
plate theory due to Mau (1973) [Journal of Applied Mechanics-Transactions of the ASME 
40, 606-607], as we treat each layer independently and enforce the compatibility by the La- 
grange multipliers. The finite clement discretization at the layer level adopts the finite-strain 

u 

shear deformable element proposed by Ibrahimbegovic and Frey (1993) [International Jour- 
nal for Numerical Methods in Engineering 36, 3239-3258]. The resulting system is solved 
by the Newton method with consistent linearization. By comparing the model predictions 
against available experimental data, analytical methods and two-dimensional finite element 
t-H simulations, we demonstrate that the proposed formulation is reliable and provides accu- 

racy comparable to the detailed two-dimensional finite element analyses. As such, it offers a 
convenient basis to incorporate more refined constitutive description of the interlayer. 
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1 Introduction 

# > 

Due to present trends in architecture and photovoltaics, the use of structural glass expands from 
traditional window panes to large-area surfaces, roof and floor systems, columns or staircases. 
This leads to an increased emphasis on the safety and the mechanical performance of structural 
members made of glass. Perhaps the most popular material system meeting these criteria is 
laminated glass. It is a composite structure produced by bonding multiple layers of glass together 
with a transparent plastic interlayer, typically made of polyvinyl butyral (PVB) foil [10\ 125] . 
The interlayer absorbs the energy impact, thereby resisting the glass penetration, and keeps the 
layers of glass bonded when fractured. 

Laminated glass units exhibit a complex mechanical response, as a consequence of their het- 
erogeneity. Namely, the contrast in elastic properties between glass and the interlayer typically 
exceeds three orders of magnitude, which renders classical laminate theories inapplicable since 
the glass layers deform mainly due to bending, whereas the interlayer experiences a pure shear. 
Moreover, PVB is a viscoelastic material exhibiting a high sensitivity to temperature changes, 
e.g. [5]. Finally, glass structures are very slender and must be analyzed using geometrically 
non-linear theories. 
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In practical applications, the behavior of laminated glass units is often approximated by an 
equivalent, geometrically linear, single-layer elastic system. According to the performance of the 
interlayer, we distinguish the layered case, in which the structure responds as an assembly of 
independent layers, and the monolithic model with thickness equal to the combined thickness 
of glass layers and interlayer. This approach was pioneered by experimental studies of Behr 
et al. [3], who demonstrated that the degree of coupling due to the interlayer is significant 
at room temperatures and ceases at elevated temperatures. This was extended later in [3] 
to quantify the validity of the monolithic approximation in terms of temperature range and 
load duration. Norville et al. [22] further refined these results by studying the shear coupling 
around the transition temperature of the interlayer and demonstrated that the performance of 
laminated units exceeds the layered limit even above the transition point. Vallabhan et al. |27j 
introduced a notion of the strength factor as a ratio between the maximal principal strength in 
the equivalent monolithic unit and the laminated system and demonstrated that it can reach 
values smaller than one for cases of practical interest. They attributed this to the effects of 
geometrical non-linearity which become significant earlier for the laminated units than for the 
monolithic ones. These developments have been recently put on a rigorous basis by Galuppi 
and Royer-Carfagni [8], who derived explicit variationally-based formulas for deflection- and 
stress-equivalent thicknesses of monolithic systems. 

Apart from experimental results, the validity of simplified models was assessed by several 
analytical studies. When restricting our attention to planar glass beams, the first study is due to 
Hooper [12] , who derived a three-layer model under small deflections and presented the solution 
for the four point bending setup with different duration of loading and ambient temperatures. 
Later, A§ik and Tescan [2] extended this model to account for large deflections and demonstrated 
that they significantly contribute to the overall response when the normal forces start to develop, 
see also Section[5]for a concrete example. In the linear setting, Ivanov [17J proposed a procedure 
for thickness optimization of triplex glasses and Foraboschi [7] demonstrated that the empirical 
rules proposed by Behr et al. [3J may lead to unsafe designs. Recently, Schultze et al. [23] 
performed an experimental-analytical study into the response of the laminated structures used 
in photovoltaic applications under three-point bending. 

Even though the analytical approaches give valuable insight into the behavior of laminated 
glass structures, they still suffer from two major limitations. First, the only closed-form solutions 
we are aware of hold only in the absence of membrane effects. In the opposite case, the governing 
differential equations have to be discretized anyhow and the discrete problem is often solved 
using the problem-specific procedures that experience convergence problems, e.g. [Tj. Second, 
the analytical approaches require the interlayer to be replaced by an equivalent elastic material, 
with properties adjusted to the loading duration and ambient temperature. As shown by Galuppi 
and Royer-Carfagni [9], this approach leads to significant errors in local stresses and strains. 

With these limitations in mind, we proposed in [28] a numerical approach to the analysis of 
laminated glass beams based on a refined laminate theory due to Mau |21j . In this framework, 
the structure is seen as an assembly of shear-deformable layers with independent kinematics and 
the inter-layer interaction is enforced by the Lagrange multipliers. Due to the modular format, 
each layer can be discretized by appropriate finite elements, and the resulting system is ideally 
suited for the deployment of efficient iterative solvers [19] or multi-scale modeling at the layer 
level [26]. 

In this paper, we extend our previous results [28], valid in small strains, to the finite-strain 
regime. To this purpose, each layer is modeled by the Reissner beam theory [23J, briefly sum- 
marized in Section [2] in a variational format. Discretization in Section [3] is based on a robust 
finite element formulation proposed by Ibrahimbegovic and Frey in [13] and later applied, e.g., 
to discrete materials modeling [13] . or optimal control and design of structures |15j. In Sec- 
tion [4j we derive the discretized system arising from the optimality conditions of the associated 
constrained optimization problem, and perform its solution by the Newton method. Accuracy 
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of the implementation is examined in Section [5] by verifying our results against data presented 
by A§ik and Tescan [2]. To make the paper self-contained, in Appendix [A| we collect the details 
on tangent operators needed when implementing the Newton method. Note that this paper is 
focused on the immediate elastic response of laminated structures. The effect of time-dependent 
viscous properties of the interlayer will be discussed independently and in more details in a 
forthcoming publication. 

The following nomenclature is used in the text. Scalar quantities are denoted by lightface 
letters, e.g. a, and the bold letters are reserved for matrices, e.g. a or A. A T standardly stands 
for the matrix transpose and A -1 for the matrix inverse. The subscript in parentheses, e.g. a®, 
is used to emphasize that the variable a is associated with the i-th. layer. 



2 Model formulation 

In our setting, a glass beam is considered as an assembly of three beams of identical length L, 
with the cross-section dimensions b x h®. Each layer is considered to behave according to the 
Reissner finite-strain beam theory [23], i.e. we assume that cross-sections remain planar, but 
not necessarily perpendicular to the deformed beam curve, and that the distance of a point at 
the cross-section from the centerline remains unchanged. For the i-th layer, the coordinates of 
a point in the deformed configuration can be determined as, Figure [TJ 

xW(I< ! U (l) ) = Of + X« +uf(X^) + sm(<p®(X®))Z®, (la) 
z®(X®, Z«) = Of + wf(X^) + cob(<p®(X®))Z®, (lb) 

(i) (i) Cil 

where O J and O z stand for the coordinates of the beam origin, X"> is the ordinate of a cross- 

(i) (i) ("\ 

section, Uq and Wq are centerline displacements measured in the global coordinate system, ipW 
is the rotation of the cross-section, and 2" is the coordinate measured along the cross-section. 
The inter-layer interaction is ensured via the geometric continuity conditions at the interfaces 
between the layers (i = 1,2) 

x®(X®, - x (m) (X», -\h^) = 0, (2a) 

zWprW, \h®) - -\h^) = 0. (2b) 



y,Y x,X 




Figure 1: Kinematics of a cross section of the bottom layer of a laminated beam (i = 3). 

As demonstrated, e.g. by Ibrahimbegovic and Frey |14j and Irschik and Gerstmayr |16j . 
the Reissner beam theory can be consistently derived from the continuum framework by the 
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Biot-type strain tensors for the kinematics parametrization Q. To this purpose, we use the 
deformation gradient in the form 

q x (i) Qx( l ) 



with individual entries provided by 

d4°(* (i) ) 



dX,® dZ® 
dz® dz® 



(X®,Z®), 



(3) 



11 ' dXW 
F«=sin(^)(xW)), 

i if/ 1 



dX« 



F«=cos(^)(X«)). 
The deformation gradient F^' admits a multiplicative decomposition 

F (i) ( X (0 , Z« ) = flW (X« ) E7« (X« , Z« ) , 
where the rotation matrix of the i-th. layer is given by, recall Figure [TJ 



cos I <p 



(c^)(X«)) sin(<pW(X®)) 



and 17 ® is follows from the inverse relation 

U®(X®,Z®) = [R {i ^{X^)Y l F (i \x^\Z^) 

that provides 

I 

dXW 



(4a) 
(4b) 

(4c) 
(4d) 

(5) 
(6) 

(7) 



IP ' 


= COS 


u {i) 

u 12 


= 0, 


u {i) 

u 21 


= sin 


u (i) 

u 22 


= 1. 



( v e> ( x«>))(i + <<- Y,,) > 



dX® 



(0 



(8b) 



(8c) 
(8d) 



Employing the definition of the Biot-type strain tensor, e.g. [18, Eq. (24.27)], 

H®(X®, Z«) = (7 (i) (X«, ZW) - J, (9) 

the non-zero strain components are given by 

H$(X®,Z®) = E^(X^) + K®(X®)Z®, H^(X^) = r®(X®), (10) 

where E^\ and denote the generalized normal strain, pseudo-curvature and shear strain 
introduced by Reissner |23| : 

E^=E{ut\wf^) 



= cos(^(X«))(l + 

r« = r(4^4V i} ) 



sm 



*0 

dX« 

.(<)/ 



sin 



K« = Jf Y(4 i) ,4 i) ,^)) 



dX« 



(11a) 

(lib) 
(11c) 
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It is useful to express these relations in a compact form 



w 

o 
o 



(X«), 



(X«), 



(12) 



and denote by -E(w) a mapping assigning the generalized strain measures to the generalized 



centerline displacements u according to Eq. (11) 



The model is completed by specifying the energy functionals associated with the deformation 
of the laminated beam. In particular, the internal energy of the i-th layer is provided by |14j 



ng (j5«>) = (is) 

pL 

i / E^W^^X^^ + G^U^rW^ 
Jo 

where and denote the Young and shear moduli, and A® = As*"* = %A^> and 

J(0 = ±b(h^) 3 stand for the the cross-section area, effective shear area, and the second moment 
of area, respectively. The external energy due to loading acting at the i-th layer assumes the 
form 



w, 



w (x«)/«(x«)dx« -E^W^W), 

P =i 



Mi) 



(14) 



where, for simplicity, we assume that the structure is subjected to the distributed loading with 

intensity and to n p ^ concentrated forces acting at Xp\ The total energy of the i-th 
layer is then given by 

n«( M ( i ))=n«(£( M »))+n»( w «), (is) 

and at the level of the whole structure it reads 



(16) 



i=l 



3 Finite element discretization 

It is convenient for the numerical treatment to discretize all layers identically with n e elements 

(i) 

per each layer; by £l e we denote the e-th element of the i-th layer. The displacement fields at 
the element level are approximated as 

u^(x) ps N e (x)df l \ w$(x) ps N e (x)dJ^, V? W (x) « N e {x)df i) for x G ft®, (17) 

where N e is the matrix of piecewise linear basis functions, and e.g. dt = [<fg\, ^e\] T collects 
the nodal rotations tp®, cf. Figure [2] 
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Figure 2: Scheme of discretization and element degrees of freedom. 
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Inserting the approximate fields (17) into the energy functionals (15) results in 



n£> («<«>) -EC(df»), !&(«<•■>) - E n ele(4°) = -£4 i)T /S2U (is) 



e=l 



e=l 



e=l 



where the matrix d e collects the nodal unknowns, 



(i) (i) (i) {€) (i) (i) 



(19) 



recall Figure 



(i) 

and /g xt e stores the corresponding generalized nodal forces. As demonstrated 



in detail by Ibrahimbegovic and Frey |14j . the normal and shear locking even for the non-linear 
kinematics can be suppressed by adopting the selective one-point integration of the corresponding 



terms in ( 13 ). Thus, the approximate generalized strain measures ( 11 ) are taken as element-wise 
constants: 



Tf{df) = * (4° + Ati®) sin/3« + 4yA<cos/3«, 



(0 



(20a) 
(20b) 

(20c) 



with, e.g., A^e = </9g*2 — <Pe\, f$ = gC'/'ei + <Pe]s)i an d denoting the element length. The 
contribution of the e-th element to the internal energy simplifies to 

n« e (d«) = i (E«A«tf«(d«) 2 + G^A^T i i\dff + E«/»K«(d«) 2 ) Lf. (21) 



The discretization is completed by enforcing the inter-layer compatibility conditions ( 2a ) at 
the element nodes, 



„(*.*+!) 



0, 



„(*»*+!) 



0, 



indexed at the level of a single layer by j = 1, 2, . . . , n e + 1, cf. Figure [TJ and with 

C g +1 ) = _J(fc(0 + + w ® - + cob V W + fc<*«) cos^ i+1) ) 

For later reference, these are introduced in a compact form as 

c(d) = 0, 



(22) 

(23a) 
(23b) 

(24) 



where d = [d^\ dP\ d^} is a 9(n e + 1) x 1 column matrix of nodal degrees of freedom and c 



collects the 4(n e + 1) compatibility conditions (22) 



4 Governing equations 

The true nodal displacements d* follow from the minimization of the discretized energy func- 
tional 



3 n c 



n(d) = EE n S, e (4 ) ) + nS t , e (4 l) ), 



(25) 



t=l e=l 
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subject to both kinematic constraints and compatibility conditions (22). While the kinematic 
constrains are dealt with by matrix reduction techniques, e.g. |18} Appendix A] , the compatibility 
is enforced via the Lagrange multipliers [21[ [26] . This yields the Lagrangian function in the form 



4(n e +l) 



C(d,X)=U{d) + \ T c(d)=U(d)+ ^mCm(d) 



(26) 



m=l 



where A is a 4(n e + 1) x 1 matrix storing the corresponding Lagrange multipliers. 

The corresponding Karush-Kuhn- Tucker optimality conditions read, e.g. [HI Chapter 14], 



V d £ (tf , A*) = VII(tf ) + Vc(d*) T A* = 0, 
V x C(d*,X*) = c{d*) = 0, 



(27a) 
(27b) 



with V denoting the gradient operator and V, designating the partial derivative with respect 
to variable •. These represent a system of non- linear equations, to be solved using the Newton 
iterative scheme. 

To that end, we assume that displacements at iteration ( k d, k X) are known and search for 
the iterative correction in the form 



fc+i 



d 



"d + k+1 8d 



The values of k+1 5d are obtained by means of the linearized expressions 

VU( k+1 d) « VU( k d) + V 2 U( k d) k+1 Sd, 
c( k+1 d) ps c{ k d) + Vc{ k d) k+1 5d, 



Vc( k+1 d) ps Vc(*d) + V 2 c( k d) k+i 5d 



2/h j\fc+l. 



(28) 

(29a) 
(29b) 
(29c) 



which, when introduced into the optimality conditions (27a) and (27b), result in a linear system 
of equations, cf. [6] and 



(30) 







' k+l 5d' 




~ k f. — f 

J int J cxt 


k C 




. fc+lA _ 




k C 



Here, we employ the short-hand notation 



4(n c +l) 



'K = V 2 U( k d)+ k ^mV 2 c m ( k d) = K t ( k d) + K x ( k d, k \), 



m=l 



k C = Vc( k d), 



VU( k d). 



(31a) 

(31b) 
(31c) 



It follows from the specific form of the energy function (25) that the stiffness matrix Kt and 
the matrices of internal / int and external forces / ext exhibit the block structure 



(32) 



For the i-th layer, these are obtained by the assembly of contributions of the e-th element 
provided by 









r>(ir 

J int 




r f (i)i 

J ext 


K t (d) = 


KfV 2) ) 


i f int 


f (2) 
J int 
f (3) 
J int_ 


) f ext 


f (2) 
J ext 
f (3) 
J ext. 



M) _ alI int,e 



(0 



u li int,e 



df (i) 

U J mt,e 



dd, 



dd] 



(33) 
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In order to keep the paper self-contained, explicit expressions for the matrices needed to set 
up the linearized system ( 30 ) are summarized in Appendix [Aj Termination of the iterative 
process (30) is driven by two residuals [6j Section 14.1] 



m 



k -£ -P i kfTk \ I 

T int ~~ Jext ~>~ ^ A \ 



cxtll2 



" c \\2 



mm 



h(*y 



(34) 



quantifying the validity of the first-order optimality conditions (27). This provides the last com- 



ponent of the non-linear iterative solver, the implementation of which is outlined in Algorithm]!] 



Algorithm 1: Conceptual implementation of the Newton method. 

Data: initial displacement °d, tolerances t\ and €2 
Result: d*, X* 

k <r- 0, °A <r- 0, assemble k f int , k c and k C 
while ( k r]i > e\) or ( k 7]2 > £2) do 
assemble k K 

solve for ( k+1 5d, k+1 X) from Eq. Q 
k+1 d^ k d + k+1 5d 
assemble fc / int , k c and k C 
_ k <— k + 1 

k d, X* <- k X 



d 



5 Examples 

In this section, the proposed finite element formulation is verified and partially validated against 



two benchmarks after A§ik and Tescan [2], involving simply supported, Section 5.1, and fixed- 



end, Section [572| three-layer beams subjected to three-point bending. The comparison is based 
on the centerline deflections and the extreme stress values at the bottom surface of the bottom 
layer (i = 3). The latter are obtained from the element-wise constant strain measures in the 
form 

S(3)= E (3)(£(3) + l^(3 )/t (3)^ (35) 



cf. Eq. (10) and Eq. (20), and extrapolated to nodes by the least-square method assuming 
the piecewise linear distribution of stresses, e.g. |11| . Accuracy of our results against reference 
experimental data is quantified by 

^exp= ( ' ) ^" ( ' ) , (36) 

and a similar approach is adopted for reference data obtained by analytical model (an) or two- 
dimensional finite element simulations (num), both from [2]. 

Results of the finite-strain formulation are complemented with the small-strain version [28J 
that corresponds to the first iteration of Algorithm [TJ For the completeness, we also provide 
the response of the equivalent monolithic beam of total thickness (h^ + + h^) and the 
layered beam corresponding to two independent layers of thicknesses and , under the 
assumption of geometric linearityj^] Note that to avoid confusion with the terminology used for 
composite structures |21j . the layered approximation will is referred to as two-layer and the 
present formulation is denoted as refined in what follows. 



1 the present section are reproducible with MATLAB® scripts available at the arXiv version of this paper. 
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5.1 Simply supported beam 



The first example concerns the simply supported beam with a span of 0.8 m and the total length 
of 1 m, with a symmetric layer setup of 5-0.38-5 mm in thickness. The structure is subject to a 
concentrated force at the mid-span ranging from 50 to 200 N and the quantities of interest are 
registered at the bottom layer, cf. Figure |3^a). The material data for individual layers appear 
in Table [T] The problem was discretized with n e = 40 elements per layer (half if symmetry is 
exploited), in order to achieve the four-digit accuracy of the mid-span deflection for F = 200 N, 
and the tolerances of the iterative solver were set to e± = 62 = 10~ 6 . 



Table 1: Material data for simply supported beam, after [2]. 
Property Value 
Young's modulus of glass, = 64.5 GPa 

Shear modulus of glass, GW = G( 3 ) 26.2 GPa 
Young's modulus of PVB, E^ 2 ) 3.61 MPa 

Shear modulus of PVB, G^ 2 ) 1.28 MPa 



The centerline deflections as predicted by the considered models appear in Figure |3J 
In this case, the deflection of the refined formulation remains bounded by the two-layer and 
monolithic cases. The response of the laminated structure is in fact closer to the monolithic beam 
than to the layered approximation, thereby demonstrating that the interlayer provides sufficient 
coupling to achieve the composite action. We also observe, in agreement with assumptions of 
many analytical models discussed in Section[TJ the absence of non-linear effects due to statically- 
determinate character of the test. 



beam width 0.1 m 



7T 

0.1 m 



0.4 m 



a) 



A 



-/- 



i— glass 5 mm 
PVB 0.38 mm 
glass 5 mm 



0.4 m 



strain gage 



7T 

0.1 m 

-X *■ 



(b) 



I'll 

^^^^ Monolithic 


Two-layer ^ Refined 
s v Linear ~ 

1,1,1," 


\ - 


* : 

non-linear 



0.1 



0.2 0.3 

Distance [ml 



0.4 



0.5 



Figure 3: Simply supported three-layer beam; (a) experiment setup and (b) centerline deflections 
for F = 50 N. 



These findings are further supported by Table [2] collecting the numerical values of mid-span 
deflections and experimental data. The results confirm that predictions of the numerical scheme 
are in a perfect agreement with the analytical model and correlate well with the experimental 
data. The limit cases, on the other hand, are too apart to be of practical use and exhibit errors 
that far exceed the difference between the refined models and experiments. 

The response of the structure to an increasing force expressed in the terms of the deflections 
and stresses appears in Tables [3] and |4j respectively. The data clearly demonstrate that the 
response of the structure remains linear even for larger loading; the differences appear only for 
F = 200 N and are indeed negligible. The deflections remain sufficiently accurate with respect 
to both the analytical model and experiments, while stress values display larger discrepancies. 
We attribute the ~ 1.3% difference between the analytical and numerical models to the fact 
that A§ik and Tescan in [2j assume glass to deform exclusively in bending and the PVB layer 
in shear only, while the present approach accounts for both effects simultaneously in all layers. 
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Table 2: Mid-span deflections Wq at = L/2 for the simply supported beam subjected to 
F = 50 N. 



Model 


Deflection [mm] 


Vexp [%] 




Experiment^ 


1.27 


X 


-5.2 


Analytical model^l 


1.34 


5.5 


X 


Laminated linear 


1.34 


5.6 


0.1 


Laminated non-linear 


1.34 


5.6 


0.1 


Monolithic 


0.89 


-30.1 


-33.8 


Layered 


3.97 


213 


196 



This error is still significantly smaller compared to the one measured against the experimental 
data. This might be attributed to inaccuracies in the experiment as explained in detail in [2]. 



Table 3: Comparison of mid-span deflections Wq at = L/2 for the simply supported beam. 



Load [N] 


Reference data^ 




Linear model 




Non-linear 


model 




w cxp [mm] 


w an [mm] 


w [mm 


] % xp [%} n. 


«, [%) 


w [mm] 


Vcxp [/ 


M Tfa [%] 


50 


1.27 


1.34 


1.34 


5.6 


0.1 


1.34 


5.6 


0.1 


100 


2.55 


2.69 


2.68 


5.2 


-0.3 


2.68 


5.1 


-0.3 


150 


4.12 


4.03 


4.02 


-2.3 


-0.1 


4.02 


-2.5 


-0.3 


200 


5.57 


5.38 


5.37 


-3.7 


-0.3 


5.35 


-4.0 


-0.6 



Table 4: Comparison of extreme stresses at X^ = L/2 for the simply supported beam. 



Load [N] 


Reference data^ 




Linear model 




Non-linear mo 


del 




Sexp [MPa] 


San [MPa] 


S [MPa] rfexp [%] 


7?an [%] 


S [MPa] 


Vexp [%] 


r/an [%] 


50 


9.55 


7.23 


7.14 


-25.3 


-1.3 


7.14 


-25.2 


-1.3 


100 


12.34 


14.45 


14.27 


15.6 


-1.3 


14.28 


15.7 


-1.2 


150 


21.89 


21.68 


21.40 


-2.2 


-1.3 


21.42 


-2.2 


-1.2 


200 


26.27 


28.90 


28.53 


8.6 


-1.3 


28.55 


8.7 


-1.2 



5.2 Fixed-end beam 

Effects of geometric non-linearity are illustrated by means of a thin 1.5 m long three-layer beam of 
thicknesses 2.12-0.76-2.12 mm subjected to a concentrated force at the mid-span with intensity 
ranging from 15 to 150 N, Figure [5~2^ a). The tolerances for the Newton method were set to the 
same value as in the previous example, and so were the material constants for the glass layers. 
For the PVB interlayer, we used E< 2 ) = 2.8 MPa and G^ 2 ) = 1 MPa. The reference data from [2] 
include the results of the analytical model, obtained by a finite-difference iterative solver [1], 
and of detailed two-dimensional large-deformation finite element simulations. To make a direct 
comparison to these values, we employ the same number of elements per layer, n e = 150, but 
analogous results are obtained for coarser discretizations. 

The centerline deflections for geometrically linear and non-linear refined beam formulations 
appear in Figure |5.2[ b) and are compared with the response of equivalent monolithic and two- 
layer formulations. While the response of the linear model still falls within the monolithic- 
two- layer bounds (with a closer proximity to the monolithic beam), the deflection of the fully 
non-linear model is considerably smaller that for the monolithic case. Moreover, the gap between 
the monolithic and two-layer cases becomes even more pronounced than in the previous example. 

The detailed numerical values presented in Table [5] further support our findings. We see that 
the layered assumption as well as small-strain hypothesis are too conservative and lead to highly 
inefficient designs, as their errors exceed 100%. The accuracy of the monolithic approximation 
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Figure 4: Fixed-end three-layer beam; (a) experiment setup and (b) centerline deflection for 
F = 15 N. 



Table 5: Mid-span deflections Wq at = L/2 for the fixed beam subjected to F = 15 N. 



Model 


Deflection [mm] 


Vnum [%] 


»7an [%) 


2D finite element model** 


5.92 


X 


0.0 


Analytical model^l 


5.92 


0.0 


X 


Laminated linear 


14.44 


144 


144 


Laminated non-linear 


6.00 


1.3 


1.3 


Monolithic 


7.85 


32.6 


32.6 


Layered 


51.48 


770 


770 



is comparable to the simply supported setup. 

Basically the same conclusions follow from the response of the structure to an increasing 
load, Tables [6] and [7j For the largest load of F = 150 N, the error of the linear laminated 
model reaches ~ 800% for the deflections and ~ 250% for the extreme stresses. The finite-strain 
formulation remains accurate in the whole range of loading, and the errors with respect to the 
detailed two-dimensional finite element model do not up-cross ~ 1% for the deflections and 
~ 2% for stresses. In fact, it slightly outperforms the analytical model [2] in terms of the stress 
values, probably due to the more refined representation of the deformation in individual layers 
as discussed above. 

As the final check of our implementation, in Table [8] we present the convergence progress for 
the load level F = 150 N. In order to investigate reliability of the method, we expose the structure 
to the full load, instead of applying it incrementally as in [2]. The results confirm significant 
degree of non-linearity in the structural response. After the first iteration, corresponding to 
the linear model, the structure is in an out-of-equilibrium state and the layer compatibility 
is violated. In the following iterations, the residuals are gradually reduced until the ninth 
iteration, after which the method exhibits quadratic convergence in the equilibrium residual 
rji and super-linear convergence in the compatibility residual r]2- This is in a full agreement 
with available results for the Newton method, e.g. [BJ Theorem 13.6], and also explains the 
convergence difficulties of the heuristic finite-difference solver reported in [2]. 

6 Conclusions 

In this paper, we have presented an efficient approach to the analysis of the laminated glass 
beams. Based on the Mau theory for layered structures, it treats each layer separately and 
enforces the inter-layer compatibility by the Lagrange multipliers. In our implementation, 
we utilize a reliable finite element formulation of the Reissner finite-strain beam theory due 
to Ibrahimbegovic and Frey to discretize individual layers, and solve the resulting system of 
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Table 6: Comparison of mid-span deflections Wq at = L/2 for the fixed-end beam. 
Load [N] Reference data^ Linear model Non-linear model 



w an [mm] w num [mm] w [mm] r? an [%] r? num [%] w [mm] r? an [%] r/ num [%] 



15 


5.92 


5.92 


14.44 


143.9 


143.9 


6.00 


1.3 


1.3 


30 


8.10 


8.10 


28.88 


256.6 


256.6 


8.17 


0.8 


0.8 


45 


9.60 


9.60 


43.32 


351.3 


351.3 


9.66 


0.6 


0.6 


60 


10.78 


10.78 


57.76 


435.9 


435.9 


10.83 


0.5 


0.5 


90 


12.64 


12.63 


86.65 


585.5 


586.0 


12.68 


0.3 


0.4 


120 


14.10 


14.09 


115.53 


719.4 


719.9 


14.14 


0.3 


0.3 


150 


15.34 


15.32 


144.41 


841.4 


842.6 


15.36 


0.1 


0.3 


Table 7: Comparison of extreme stresses at A^ 3 ) 


= L/2 for the fixed-end beam. 




Load [N] 


Reference data^ 




Linear model 


Non-linear model 




San [MPa] 


5 num [MPa] 


S [MPa 


] ??an[%] 


??num [%] 


S [MPa] 


r?an [%] 


Vnum [%] 


15 


12.87 


12.46 


19.50 


51.5 


56.5 


12.60 


-2.1 


1.1 


30 


20.69 


19.89 


38.98 


88.4 


96.0 


20.12 


-2.7 


1.2 


45 


27.13 


25.94 


58.44 


115.4 


125.3 


26.28 


-3.2 


1.3 


60 


32.82 


31.25 


77.88 


137.3 


149.2 


31.69 


-3.4 


1.4 


90 


42.82 


40.51 


116.69 


172.5 


188.1 


41.18 


-3.8 


1.6 


120 


51.68 


48.64 


155.43 


200.8 


219.6 


49.53 


-4.2 


1.8 


150 


59.76 


56.00 


194.10 


224.8 


246.6 


57.13 


-4.4 


2.0 



equations iteratively by the Newton method with consistent linearization. On the basis of the 
performed simulations, we conjecture that 

• in the absence of membrane effects, the formulation reduces exactly to the small-strain 
model introduced in our previous work, 

• although the discretization is based on the lowest-order polynomial basis functions, the 
method provides results with accuracy comparable to the detailed two-dimensional large- 
strain finite element simulations, 

• the Newton method exhibits a reliable super-linear convergence even for high degrees of 
non-linearity. 

Extension of the current framework to include temperature- and time-dependent properties 
of the interlayer will be reported independently. 
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A Sensitivity analysis 



The expression for the internal nodal forces follows directly from Eq. (33)i. After certain ma- 
nipulations we arrive at, cf. |20| 
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By an analogous procedure one obtains from Eq. (33)2 
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The remaining terms in the system ( 30 ) originate from the compatibility conditions ( |22[ ) . In 
particular, the matrix C is analogous to the small-strain tying condition [281 Section 4]. The 
block of C, associated with a node j and layers i and (i + 1) attains the form 
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The second derivatives of the compatibility conditions quantify their contributions to the tangent 
stiffness 
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This additional term is expressed as 
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with non-zero entries provided by 
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